3. Mathematical Background

In order to solve the obtained PDEs (2.5) or (2.9) various numerical methods are developed in this toolbox. Users can simulate the models with these methods to evaluate the results.

Actually, first, (2.5) and (2.9) are discretized in time by backward-Euler method and the following equations are obtained [MRP21]:

(3.1)\[\frac{\rho(\theta,t_{i})-\rho(\theta,t_{i-1})}{\Delta t}=-\frac{\partial}{\partial\theta}\Big((\omega+ u(t)^T\mathcal{Z}(\theta))\rho(\theta,t_{i})\Big)~~~~i=0, ..., m ,\]

(3.2)\[\frac{\rho(\theta,t_{i})-\rho(\theta,t_{i-1})}{\Delta t}=-\frac{\partial}{\partial\theta}\Big((\omega+ u(t)^T\mathcal{Z}(\theta))\rho(\theta,t_{i})\Big)+\mathcal{B}\frac{\partial ^2\rho(\theta,t_{i})}{\partial \theta ^2}~~~~i=0, ..., m,\]

in which, \(\Delta t\) is the time step size, \(t_i=i\Delta t\) and \(T=m\Delta t\). This scheme preserves the first-order accuracy in time. For ease of explanation, we show \(\rho^i=\rho(\theta,t_{i})\). So, (2.5) and (2.9) convert to:

(3.3)\[\rho^i+\Delta t\Big(\omega\rho^i_{\theta} +u^{i-1}(\mathcal{Z}(\theta)\rho^i_{\theta}+\frac{d\mathcal{Z}(\theta)}{d\theta}\rho^i)\Big)=\rho^{i-1}.\]

(3.4)\[\rho^i+\Delta t\Big(\omega\rho^i_{\theta} +u^{i-1}(\mathcal{Z}(\theta)\rho^i_{\theta}+\frac{d\mathcal{Z}(\theta)}{d\theta}\rho^i)-\mathcal{B}\rho^i_{\theta\theta}\Big)=\rho^{i-1}.\]

In the following sections, the proposed numerical methods are examined.

3.1. Finite difference method

The finite difference (FD) method is one of the most famous and classical numerical approaches for solving ordinary/partial differential equations. There exist several kinds of FD algorithms. We utilize the backward difference in temporal dimension and for spatial one, we apply the central and 5-point stencil to discretize PDEs (2.5) and (2.9), respectively.

Consider the following equations [For88]:

(3.5)\[ \frac{\partial\rho(\theta,t)}{\partial t}=\frac{\partial\rho(\theta,t_{i+1})-\partial\rho(\theta,t_i)}{\Delta t},~~~i=0,...,M,\]

(3.6)\[ \frac{\partial\rho(\theta,t)}{\partial\theta}=\frac{\rho(\theta_{j+1},t)-\rho(\theta_{j-1},t)}{2\Delta x},~~j=0,...,N,\]

where \(\Delta t\) and \(\Delta x\) denote the time and space steps. By imposing these formulas in (2.5), we have:

(3.7)\[ \Big(\frac{-\omega\Delta t-u^{i}\mathcal{Z}(\theta_j)\Delta t}{2\Delta x}\Big)\rho^{i+1}_{j-1}+\Big(1+u^i\mathcal{Z}'(\theta_j)\Big)\rho^{i+1}_j+\Big(\frac{\omega\Delta t+u^{i}\mathcal{Z}(\theta_j)\Delta t}{2\Delta x}\Big)\rho^{i+1}_{j+1}=\rho_j^i,\]

in which, \(i\) and \(j\) represent time and space nodes, respectively.

The obtained PDEs have periodic boundary conditions; so, the first and last points are the same (See Fig. 3.1).

Fig. 3.1 Space nodes according to the periodic boundary condition.

If we write (3.7) for the first, last, and entire points, at each time step a linear system of algebraic equations is constructed that can be solved easily with a proper solver like QR algorithm.

The derivatives in (2.9) are computed by the following formulas [Sau12]:

(3.8)\[ \frac{\partial\rho(\theta,t)}{\partial t}=\frac{\partial\rho(\theta,t_{i+1})-\partial\rho(\theta,t_i)}{\Delta t},~~~i=0,...,M,\]

(3.9)\[ \frac{\partial\rho(\theta,t)}{\partial\theta}=\frac{\rho(\theta_{j-2},t)-8\rho(\theta_{j-1},t)+8\rho(\theta_{j+1},t)-\rho(\theta_{j+2},t)}{12\Delta x},~~j=0,...,N,\]

(3.10)\[ \frac{\partial^2\rho(\theta,t)}{\partial\theta^2}=\frac{-\rho(\theta_{j-2},t)+16\rho(\theta_{j-1},t)-30\rho(\theta_{j},t)+16\rho(\theta_{j+1},t)-\rho(\theta_{j+2},t)}{12\Delta x^2},~~j=0,...,N,\]

Now, by using these equations in (2.9), we have:

(3.11)\[\frac{\mathcal{B}\Delta t}{12\Delta x^2}\rho_{j-2}^{i+1}-\Delta t\Big(\frac{\omega+u^i\mathcal{Z}(\theta_j)}{2\Delta x}+\frac{16\mathcal{B}}{12\Delta x^2}\Big)\rho_{j-1}^{i+1}+\Big(1+\Delta t u^i\mathcal{Z}'(\theta_j)+\frac{30\mathcal{B}\Delta t}{12\Delta x^2}\Big)\rho_{j}^{i+1}\\+\Delta t\Big(\frac{\omega+u^i\mathcal{Z}(\theta_j)}{2\Delta x}-\frac{16\mathcal{B}}{12\Delta x^2}\Big)\rho_{j+1}^{i+1}+\Big(\frac{\mathcal{B}\Delta t}{12\Delta x^2}\Big)\rho_{j+2}^{i+1}=\rho^i_j.\]

The obtained equations can be re-written in matrix form and solved by QR method.

3.2. Spectral method

High order numerical schemes have been used to simulate and understand the dynamics of biological phenomena such as [MD20][KAK+20][PMLR19]. In fact, spectral methods are one of the most useful tools in numerical simulations due to their exponential convergence and benefits in implementation such as the possibility of parallelization. In the pseudo-spectral collocation method, the solution is approximated by a truncated series of polynomial interpolations based on appropriate nodes. These nodes have a significant effect on the accuracy and convergence; so it is shown that an accurate approximation with exponentially convergence for the smooth problems can be achieved by using Gaussian quadrature collocation points [GPH+10].

In this toolbox, generalized Lagrange Jacobi Gauss–Lobatto functions (GLJGLFs) are considered as the basis functions. These functions have some benefits in implementation and computing complexity. We apply the generalized Lagrange Jacobi Gauss–Lobatto collocation method on (3.3) and [time_disc2_noise]. The standard Jacobi polynomials which is denoted by \(J^{(\alpha,\beta)}_n(x)\), \(n=0, 1, ...\) are one of the most useful orthogonal polynomials on [-1,1]. Besides, these polynomials can be constructed via following recursive formula [HRBPR19]:

(3.12)\[ J_0^{(\alpha ,\beta)}(x)=1,~~~~~J_1^{(\alpha ,\beta)}(x)=\frac{\alpha +\beta+2}{2}x+\frac{\alpha -\beta}{2},\\ J^{(\alpha ,\beta)}_{n+1}(x)=(a^{\alpha,\beta}_{n}x-b^{\alpha,\beta}_{n})J_{n}^{(\alpha ,\beta)}(x)-c^{\alpha,\beta}_{n}J_{n 1}^{(\alpha ,\beta)}(x)~~~n\geq 1\]

where

(3.13)\[a^{\alpha,\beta}_{n}=\frac{(2n+\alpha+\beta+1)(2n+\alpha+\beta+2)}{2(n+1)(n+\alpha+\beta+1)},\\ b^{\alpha,\beta}_{n}=\frac{(\beta ^2-\alpha ^2)(2n+\alpha+\beta+1)}{2(n+1)(n+\alpha+\beta+1)(2n+\alpha+\beta)},\\ c^{\alpha,\beta}_{n}=\frac{(n+\alpha)(n+\beta)(2n+\alpha+\beta+2)}{(n+1)(n+\alpha+\beta+1)(2n+\alpha+\beta)}.\]

The Jacobi polynomial \(J_n^{\alpha,\beta}(x)\) has exactly n real roots on [-1,1].

Generalized Lagrange Jacobi functions are an extension of classical Lagrange polynomials based on Jacobi Gaussian nodes.

We define GLFs corresponding to the interpolation nodes \(x_0<x_1<...<x_n\) as follows:

(3.14)\[L_i^{\phi}(x)=\prod_{j=0}^n \frac{\phi(x)-\phi(x_j)}{\phi(x_i)-\phi(x_j)},\]

in which \(\phi(x)\) is an arbitrary continuous and sufficiently differentiable function depending on the problem features. According to the function \(\phi(x)\) different basis functions can be defined in the various intervals.

Due to the properties of GLJGLFs, their derivatives can be computed by a direct formula; therefore, there is no need to derive directly from functions, which reduces the computational difficulty and cost dramatically.

We can state two main advantages for this method:

• First, \(\phi(x)\) is an arbitrary function that can be selected according to the domain and physical properties of the problem. A clever choice increases the effectiveness and accuracy. With different choices of \(\phi(x)\), many new basis functions can be generated at different intervals, for instance [DP19]:

  • Classical Lagrange functions: \(\phi(x)=x\).

  • Shifted Fractional Lagrange functions on [0,L]: \(\phi(x)=2(\frac{x}{L})^{\alpha}-1\).

  • Rational Lagrange functions on \([0,\infty): \phi(x)=\frac{x-L}{x+L}\), where L is a positive real constant.

  • Exponential Lagrange functions on \((-\infty,\infty):\phi(x)=e^x\).

• Since derivative matrices can be obtained by direct formulas, there is no need to calculate the derivatives of the functions. Moreover, generalized Lagrange functions have the property of the Kronecker delta; therefore, the proposed method has a small computational cost.

In our problems, we need first and second derivatives matrices and we denote them as \(\mathbf{D}^{(1)}\) and \(\mathbf{D}^{(2)}\), respectively. These matrices are pre-calculated and stored in the matrices folder. For more information on how to calculate them, interesting readers can refer to [PLDM18][PLMD18][LD20][MHRL+20].

Since in our model \(\theta \in [0,2\pi]\), we need to change the domain of Lagrange polynomials; thus, we consider \(\phi(x)=\frac{x}{\pi}-1\).

In the pseudo-spectral method, one can approximate function \(\rho^i\) for time step i in (3.3) and (3.4) as follows [MRP21]:

(3.15)\[\rho^i(\theta)\approx \hat{\rho}^i(\theta)= \sum_{k=0}^n\rho^i(\theta_k)L_k^{\phi}(\theta)=\mathbf{L}^{\phi}(\theta)\mathbf{P}^i\]

where

(3.16)\[\mathbf{L}^{\phi}(\theta)=\begin{bmatrix} L_0^{\phi}(\theta),~L_1^{\phi}(\theta),~ ... ,~ L_n^{\phi}(\theta) \end{bmatrix},\\ \mathbf{P}^{i}=\begin{bmatrix} \rho^i(\theta_0),~\rho^i(\theta_1),~ ... ,~ \rho^i(\theta_n) \end{bmatrix}^T,\]

The Gauss-Lobatto-Jacobi nodes, \(\{\theta_k\}_{k=0}^n\), are chosen as the collocation points to reduce the computations and satisfy the boundary conditions. Now, by using the aforementioned matrices, the problems (3.3) and (3.4) are reduced to linear systems of algebraic equations which can be solved with direct solvers as QR or LU factorization.

3.3. Radial basis function method

Meshfree methods are another category of methods for solving differential equations. Compared to other numerical methods, meshfree methods have a low computational cost because they do not require a pre-designed mesh in the problem domain [SVLH18][DA19][AD19][JS18][SJ18b][SJ18a]. Indeed, in these methods, a set of scattered data represent the domain and boundaries of the problem [Fas07]. One of the most popular meshfree methods is the method based on radial basis functions (RBFs) that have made remarkable progress over the last decade. High execution speed, ease of implementation, high accuracy, and flexibility are prominent features of this method [Fas05][Sar05][SL16][SVHL15][SMN+15][RHookLS18].

In this toolbox, two meshfree methods have been used to numerically simulate the control model (3.3) and (3.4). One is the collocation method based on RBFs and the other is the radial basis function generated finite difference method (RBF-FD) which will be explained in the next section. In the RBF collocation method, various radial basis functions can be used, which we have used the Wendland compact support function [Wen95] here:

(3.17)\[\phi_i(x)=(1-\delta r_i)^8_+(1+8\delta r_i+25\delta r_i^2+32 \delta r_i^3),\]

where \(r_i=\|x-x_i\|\) is the distance from node \(x_i\) to \(x\), while \(\delta\) is the size of support for radial function \(\phi_i(x)\). In addition, the sign ‘+’ in the function means that \((1-\delta r_i)^8_+\) is \((1-\delta r_i)^8\) for \(\delta r_i \in [0, 1)\) and vanished in other regions.

The following definitions are the main features of RBF functions.

[Fas07] A function \(\phi : \mathbb{R}^s \rightarrow \mathbb{R}\) is called radial provided there exists a uni-variate function \(\varphi [0, \infty)\rightarrow \mathbb{R}\) such that \(\phi(\textbf{x})=\varphi(r),~ where=\|\textbf{x}\|\) and \(\|.\|\) is some norm on \(\mathbb{R}^s\), usually the Euclidean norm.

In RBF collocation method, $rho^i$ for time step i in (3.3) and (3.4) are approximated by linear combination of RBFs and unknown coefficients as follows:

(3.20)\[\rho^i(\theta)\approx \hat{\rho}^i(\theta) = \sum_{k=0}^n a^i_k \phi_k(\theta)=\Phi(\theta)\textbf{A}^i\]

where

(3.20)\[\textbf{A}^i=[a^i_0, a^i_1, \cdots, a^i_n], \Phi(\theta)=[\phi_0(\theta), \phi_1(\theta), \cdots, \phi_n(\theta)].\]

Derivatives of \(\rho^i\) are also easily interpolated by calculating the RBF derivative as follows:

(3.20)\[\frac{d\rho^i(\theta)}{d\theta}\approx \frac{d\hat{\rho}^i(\theta)}{d\theta} = \sum_{k=0}^n a^i_k \frac{d\phi_k(\theta)}{d\theta} \frac{d^2\rho^i(\theta)}{d\theta^2}\approx \frac{d^2\hat{\rho}^i(\theta)}{d\theta^2} = \sum_{k=0}^n a^i_k \frac{d^2\phi_k(\theta)}{d\theta^2}.\]

In this method, collocation points \(\{\theta_k\}_{k=0}^n\) can be defined in any way, as long as they cover the domain and boundaries of the equation. Specifically, in this toolbox we consider \(n+1\) nodes \(\theta_0, \theta_1, \cdots, \theta_n\) where two points \(\theta_0\) and \(\theta_n\) are located on the boundaries to apply the periodic boundary conditions, and we define \(\theta_i=(i-1)\Delta\theta, i=1,\cdots,n\) where \(\Delta\theta=\frac{2\pi}{n-1}\).

Now by placing the approximations of the function \(\rho^i(\theta)\) and its derivatives in (3.3) and (3.4), the problem becomes a set of algebraic equations in which the unknown coefficients \(\textbf{A}^i\) are obtained by solving it.

3.4. Radial basis function generated finite difference method

The RBF-FD method is the second Meshfree method used in this toolbox. This method can be considered as a generalization of the finite difference method [Tol00]. The main feature of the RBF-FD method is that it has the sparsity property of the finite-difference method and at the same time its accuracy is close to the RBF method. In RBF-FD, the derivatives of a function are approximated by the linear combination of the RBF functions at the scattered points, its framework operates like a finite difference method, so the term ”FD” is annexed to this method. By inheriting the property of FD, the coefficient matrices obtained are well-conditioned and their bandwidth becomes smaller when the number of stencil nodes decreases [FFBB16][FLP13][BMCK10].

Let us consider a set of \(n+1\) scattered points in domain \(\Omega, \Omega^s\), i.e \(\textbf{x}=\{x_0, x_1, \cdots, x_n\}\), we define a sub-domain which is the local domain associated with the point \(x_s \in \textbf{x}\). Now consider the set of indexes related to the sub-domain members:

(3.21)\[\mathcal{I}(x_s)=\mathcal{I}_s=\{j; \|x_j-x_s\|\geq r_s\}.\]

Therefore, the derivative operator \(\mathcal{L}\) in this method is calculated as follows:

(3.22)\[\mathcal{L}u(x_s)\approx \sum_{k\in \mathcal{I}(x_s)}w_i^s u(x_i),\]

where \(\{w_i^s\}\) are the unknown weights that have to be determined. It should be noted that we need to obtain the unknown vector \(\textbf{x}=\{a_i\}_{i\in \mathcal{I}(x)}\) such that \(\sum_{i\in \mathcal{I}(x)}a_i\phi(\|x-x_i\|_2)\) interpolate \(u(x)\) at each of points in the stencil \(\Omega_s\), for this purpose, we use the local distance matrix \(\Phi_s\) in the form

(3.23)\[\textbf{a}=\Phi_s^{-1}\textbf{u}_s,\]

where \(\textbf{u}_s\) is the data corresponding to stencil of \(x_s\). Likewise, the approximation of the operator \(\mathcal{L}\) is as follows:

(3.24)\[ \mathcal{L}u(x_s)= \sum_{k\in \mathcal{I}(x_s)}a_i\mathcal{L}\phi(\|x_s-x_i\|_2)=\mathcal{L}_\phi\textbf{a}=\mathcal{L}_\phi\Phi_s^{-1}\textbf{u}_s,\]

in which \(\mathcal{L}_\phi\) represents an \(1\times (n+1)\) vector containing elements \(\{\phi(\|x_s-x_i\|_2)\}_{i\in \mathcal{I}(x)}\) and \(\phi\) is a RBF centered at \(x_s\). In this toolbox, we decide to use the Gaussian RBF [BMK12] with \(C^\infty\) smoothness degrees as follows:

(3.25)\[ \phi_i(x)=e^{(-c^2r^2_i)},\]

where \(r_i=\|x-x_i\|\) is the distance from node \(x_i\) to \(x\), and the shape parameter \(c\) governs the flatness of the RBFs.

Therefore, based on (3.22) and (3.24), the following \((n+1)\times(n+1)\) linear equation is obtained:

(3.26)\[\textbf{w}^s\Phi_s=\mathcal{L}_\phi,\]

where \(\textbf{w}^s\) is a \(1\times (n+1)\) local differential weights vector corresponding to stencil \(\Omega_S\) of \(x_S\). Now, the weights of each stencil center \(x_s\); where \(s\in\{1,\cdots,n+1\}\) is obtained by forming a \((n+1)\times(n+1)\) linear system and solving it. Also by placing the weights obtained in the \(n+1\) rows, a sparse differential matrix \(\textbf{W}\) with n non-zero elements per row is obtained as follows:

(3.27)\[\textbf{W}=\begin{pmatrix} \textbf{w}^1\\ \textbf{w}^2\\ \vdots\\ \textbf{w}^{(n+1)} \end{pmatrix}.\]

Now, using weight matrices \(\textbf{W}\) and \(\textbf{W}^{(2)}\), similar to the finite difference method, the equations (3.3) and (3.4) are transformed into a set of algebraic equations, which can be solved by methods such as LU and QR, and an approximation of the function \(\rho^i(\theta)\) in a set of nodes \(\theta=\{\theta_1,\theta_2,\cdots, \theta_{n+1}\}\) at each time steps \(i\) are obtained.

3.5. Fourier decomposition method

In this approach, the probability distribution and PRC functions are approximate by finite Fourier series [MM19b]:

(3.28)\[\rho(\theta,t)\approx\frac{1}{2\pi}+\Sigma_{k=1}^{N-1}\Big(A_k\cos(k\theta)+B_k\sin(k\theta)\Big),\]

(3.29)\[\rho_f(\theta,t)\approx\frac{1}{2\pi}+\Sigma_{k=1}^{N-1}\Big(\hat{A}_k\cos(k\theta)+\hat{A}_k\sin(k\theta)\Big),\]

(3.30)\[ \mathcal{Z}(\theta)\approx C_0+\Sigma_{k=1}^{N-1}\Big(C_k\cos(k\theta)+\hat{C}_k\sin(k\theta)\Big).\]

It is worth mentioning that this method guarantees that the phase distribution is \(2\pi\)-periodic. By multiplying \(\cos(k\theta)\) and \(\sin(k\theta)\) on both sides of (3.28) and integrating from 0 to \(2\pi\) with respect to \(\theta\), we have:

(3.31)\[ A_k=\pi^{-1}\int_0^{2\pi}\rho(\theta,t)\cos(k\theta)d\theta,\]

(3.32)\[ B_k=\pi^{-1}\int_0^{2\pi}\rho(\theta,t)\sin(k\theta)d\theta.\]

Now, by taking derivative with respect to time, using integrating by parts and imposing periodic boundary condition according to the models with or without noise, we obtain:

(3.33)\[ \dot{A}_k=-k\omega B_k-\mathcal{I}_k^A(t)u(t),\]

(3.35)\[ \dot{B}_k=k\omega A_k+\mathcal{I}_k^B(t)u(t),\]

or

(3.35)\[ \dot{A}_k=-k\omega B_k-\mathcal{I}_k^A(t)u(t)-\frac{\mathcal{B}^2}{2}k^2A_k(t),\]

(3.36)\[\dot{B}_k=k\omega A_k+\mathcal{I}_k^B(t)u(t)-\frac{\mathcal{B}^2}{2}k^2B_k(t),\]

where

(3.37)\[\mathcal{I}_k^A(t)=\frac{k}{\pi}\int_0^{2\pi}\mathcal{Z}(\theta)\rho(\theta,t)\sin(k\theta)d\theta, \mathcal{I}_k^B(t)=\frac{k}{\pi}\int_0^{2\pi}\mathcal{Z}(\theta)\rho(\theta,t)\cos(k\theta)d\theta.\]

Similarly, we have:

(3.38)\[\dot{\hat{A}}_k=-k\omega\hat{B}_k,\]

(3.39)\[\dot{\hat{B}}_k=k\omega \hat{A}_k.\]

In order to approximate the functions, we should use a suitable iterative method like the Runge-Kutta algorithm to obtain the Fourier coefficients.

3.6. Runge-Kutta algorithm

At each time step in the simulation, we need to calculate the phase of the oscillators by solving (2.4) or (2.8). In this toolbox, we utilize the well-known Runge-Kutta approach which is a good choice for solving ODEs. In the following, we describe this algorithm.

3.6.1. Fourth-order Runge-Kutta algorithm for ODEs

In order to find the phase of noise free oscillators, we solve (2.4). First, we define:

(3.40)\[f(\theta)=\omega +u\mathcal{Z}(\theta),\]

then, (2.4) is discretized to the following equation (where \(k = 0,··· ,r\), and \(r\) is the number of Runge–Kutta steps, and \(h=\frac{\Delta t}{r}\)):

(3.41)\[\hat{\theta}^{k+1}_j=\hat{\theta}^{k}_j+\frac{h}{6}\bigg(\mathcal{K}_1+2\mathcal{K}_2+2\mathcal{K}_3+\mathcal{K}_4\bigg),\]

where

(3.42)\[\mathcal{K}_1=f(\hat{\theta}^{k}_j),~\sigma_1 =\hat{\theta}^{k}_j+0.5h\mathcal{K}_1,\\ \mathcal{K}_2=f(\sigma_1),~\sigma_2 =\hat{\theta}^{k}_j+0.5h\mathcal{K}_2,\\ \mathcal{K}_3=f(\sigma_2),~\sigma_3 =\hat{\theta}^{k}_j+h\mathcal{K}_3,\\ \mathcal{K}_4=f(\sigma_3).\]

where the first and last steps are \(\hat{\theta}^0_j=\theta_j^i\) and \(\theta_j^{i+1}=\hat{\theta}^{r}_j\).

3.6.2. Fourth-order Runge-Kutta algorithm for stochastic ODEs

As we know, (2.8) should be solved to find the phase of noisy oscillators at each step. We utilize the well-known fourth-order Runge-Kutta approach to solve this stochastic ODE. We have:

(3.43)\[d\theta_j^i=\delta dt+\mathcal{Z}(\theta)\bigg[u^{i-1}+\sqrt{2D}dW_j(t)\bigg],~~~j=1, \cdots, M,\]

in which \(W_j(t)\) is the standard Weiner process. We define functions \(f\) and \(g\) from (3.43):

(3.44)\[f(\theta)=\omega +u\mathcal{Z}(\theta),\]

(3.45)\[g(\theta)=\sqrt{2D}+\mathcal{Z}(\theta).\]

The main problem is discretized as follows (\(k=0, \cdots, r\), where \(r\) is the number of Runge-Kutta steps and \(h=\frac{\Delta t}{r}\)).

(3.46)\[\hat{\theta}^{k+1}_j=\hat{\theta}^{k}_j+\frac{1}{6}\bigg(h(\mathcal{K}_1+2\mathcal{K}_2+2\mathcal{K}_3+\mathcal{K}_4)+W_j(\mathcal{S}_1+2\mathcal{S}_2+2\mathcal{S}_3+\mathcal{S}_4)\bigg),\]

where

(3.47)\[\mathcal{K}_1=f(\hat{\theta}^{k}_j),~\mathcal{S}_1=g(\hat{\theta}^{k}_j),\\ \sigma_1 =\hat{\theta}^{k}_j+0.5h\mathcal{K}_1+0.5W_j\mathcal{S}_1,\\ \mathcal{K}_2=f(\sigma_1),~\mathcal{S}_2=g(\sigma_1),\\ \sigma_2 =\hat{\theta}^{k}_j+0.5h\mathcal{K}_2+0.5W_j\mathcal{S}_2,\\ \mathcal{K}_3=f(\sigma_2),~\mathcal{S}_3=g(\sigma_2),\\ \sigma_3 =\hat{\theta}^{k}_j+h\mathcal{K}_3+W_j\mathcal{S}_3,\\ \mathcal{K}_4=f(\sigma_3),~\mathcal{S}_4=g(\sigma_3).\]

Note that \(\hat{\theta}^0_j=\theta_j^i\) and \(\theta_j^{i+1}=\hat{\theta}^{r}_j\).